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Abstract 

The goal of this paper is to give confidence regions for the excursion set of a spatial function above a given threshold 
from repeated noisy observations on a fine grid of fixed locations. Given an asymptotically Gaussian estimator of 
the target function, a pair of data-dependent nested excursion sets are constructed that are sub- and super-sets of 
the true excursion set, respectively, with a desired confidence. Asymptotic coverage probabilities are determined via 
a multiplier bootstrap method, not requiring Gaussianity of the original data nor stationarity or smoothness of the 
limiting Gaussian field. The method is used to determine regions in North America where the mean summer and 
winter temperatures are expected to increase by mid 21st century by more than 2 degrees Celsius. 
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1 Introduction 


Our motivation comes from the following problem. Faced with a global change in temperature over the 
globe within the next century, it is important to assess which geographical regions are particularly at risk 
of extreme temperature change. The data used here, obtained from the North American Regional Climate 
Change Assessment Program (NARCCAP) project (Mearns et al. 2009, 2012| 20131, consists of two sets of 29 
spatially registered arrays of mean seasonal temperatures for summer (June-August) and winter (December- 
February) evaluated at a fine grid of fixed locations 0.5 degrees in geographic longitude and latitude apart 
over North America over two time periods: late 20th century (1971-1999) and mid 21st century (2041- 
2069). Specifically, the data was produced by the WRFG climate model ( |Michalakes et ah 2004) using 
boundary conditions from the CGCM3 global model (Flato 2005). We would like to determine the regions 


whose difference in mean summer or winter temperature between the two periods is greater than the 2°C 
benchmark (Rogelj et al. 2009 Anderson and Bows 20111. However, the observed differences may be 


confounded by the natural year-to-year temperature variability. Can we set confidence bounds on such 
regions that reflect the year-to-year variability in the data? 

Unlike the usual data setup of spatial statistics, the above data setup is more similar to that of population 
studies in brain imaging, where a difference map between two conditions is estimated from repeated co-located 


image observations at a fine spatial grid under those conditions (see e.g. Worsley et al. (1996), Genovese 
et al. (2002 1 , Taylor and Worsley (2007 1 , Schwartzman et al. (2010)). The methods in this paper are inspired 
by that kind of analysis. 

In general, suppose that we observe n random fields T^s), * = 1,..., n, over a spatial domain S, modeled 
as realizations of a general linear model indexed by s G S. The target function [i : S — > R. could be one of 
the parameters in the model indexed by s, in our case the mean difference temperature field. With a proper 
design, fitting the linear model at each location s will produce a consistent and asymptotically Gaussian 
estimator fi n : S —> R. as n increases. Asymptotically Gaussian estimators indexed by s also appear in 
nonparametric density estimation and regression. In those settings n would be the number of sample points. 
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Let A c be the excursion set of /i above a fixed threshold c, defined as A c := A c ( /z) := {s £ S : fi( s) > c}, 
and denote the analog for p, n by A c := A c (fi n ). We wish to obtain confidence regions Af that are nested in 
the sense that At C A c C A~ and for which the probability that 

At CA C C A~ (1) 


holds is asymptotically above a desired level, say 90%. The sets Af here are obtained as excursion sets of 
the standardized observed field /t and we call them Coverage Probability Excursion (CoPE) sets. Assuming 
that the estimated field fi satisfies a central limit theorem (CLT), we show that the probability that ([!]) holds 
is given asymptotically by the distribution of the supremum of the limiting Gaussian random field on the 
boundary dA c of the true excursion set. Using a plug-in estimate for the unknown boundary, we propose a 
simple and efficient multiplier bootstrap procedure (Wu[ 1986 Hardle and Mammen 1993 Mammen 1992 


1993), that does not require estimating the unknown (not necessarily stationary) correlation function of the 
limiting field. The validity of this procedure for very high-dimensional data has recently been shown by 


Chernozhukov et al. (2013). 


For illustration, Figure ljshows CoPE sets for the temperature data. The regions within the red boundary 
(At) have the highest confidence of being at risk, while the regions outside the green boundary (A~) have 
the highest confidence of not being at risk. Over repeated sampling, there is a probability of about 90% that 
the regions at risk include those within the red boundary and exclude those outside the green boundary. 

The problem of finding confidence sets for spatial excursion sets, sometimes also called exceedance regions 
or level sets, has been studied in the past in two major contexts that substantially differ from the problem 
under consideration here. In the geostatistics literature, the target function is itself a Gaussian field. In 
consequence, the excursion and the contour sets are random themselves. The data in this setting is a partial 
realization of the field, that is, the values of a realization of the field at relatively few spatial locations. This 
severe limitation of available information is compensated by assuming that the covariance structure of the 
field is known. This problem has been addressed from a frequentist perspective in terms of confidence regions 


for level contours (Lindgren and Rychlik 1995 Wameling 2003 French 20141 and for excursion sets (French 


and Sain 20131. Incidentally, our techniques share some similarities with French (20141, although we will 


show that distinguishing between level contours and excursion sets is important. In a Bayesian setting for 
latent Gaussian models, Bolin and Lindgren address uncertainty in both, contours and excursion sets. 

The second setting in which the problem has received attention is non-parametric density estimation 
and regression. Here, the target function is a probability density or regression function, estimated from 
realizations of a random variable with values in R 9 for some q. While the estimation of both level sets and 


contours have been well studied (Tsybakov 


2007 Singh et al. 


2009 Rigollet and Vert 


1997 Cavalier 1997 Cuevas et al. 2006 Willett and Nowak 


2009), there is less literature on confidence statements. Mason 


and Polonik (20091 showed asymptotic normality of plug-in level set estimates with respect to the measure of 


symmetric set difference. Mammen and Polonik (2013) proposed a bootstrapping scheme to obtain confidence 
sets analogous to our CoPE sets from vector-valued samples. 

The problem of finding the threshold for our CoPE sets involves computation of the tail probability of 
the supremum of a limiting Gaussian random field. In French (2014) this computation was done by Monte 
Carlo simulation assuming that the covariance structure of the field is known. More generally for unknown 
covariance function, as we attempt here, this problem was solved elegantly by Taylor and Worsley (2007) 
using the Gaussian kinematic formula. However, this method requires that the observations themselves 
be Gaussian and requires the field to be differentiable. The multiplier bootstrap allows us to avoid both 
these assumptions while being extremely fast to compute. We compare the finite sample performance of the 
Gaussian kinematic formula method and the multiplier bootstrap in a simulation. 

All computations in this paper were performed using R (R Core Team 2014). All required functions 


for computation and visualization of CoPE sets and in particular an implementation of the Algorithm [T] are 
available in the R-package cope (Sommerfeld 2015). 


Outline of the paper 

In Section [ 2 ] we propose a thresholding scheme to obtain CoPE sets Af as in (|lj) from an estimator jl of /z, 
only requiring continuity of /z and, most importantly, that fi is asymptotically Gaussian. We show that the 

















































































































2 Error control for excursion sets - CoPE sets 


3 



Fig. 1: Output of our method for the increase £>i(s) of the mean summer (June-August) and winter 
(December-February) temperatures. Shown are heat maps of the estimator &i(s). The uncertainty in 
the excursion set estimate A c = A c (bi) (purple boundary) for c = 2°C is captured by the CoPE sets 
Af (red boundary) and A~ (green boundary). The threshold was obtained according to Theorem [l] 
to guarantee inclusion Af C A c C A~ with confidence 1 — a = 0.9. The horizontal and vertical axes 
are indexed in degrees longitude and latitude, respectively. 


asymptotic coverage probability is equal to the tail probability of the limiting Gaussian field on the boundary 
dA c of the excursion set A c . 

Section [3] is devoted to presenting results and algorithms for the construction of CoPE sets when the 
target function is the parameter function in a general linear model. First, in Section |3.1| we derive central 
limit theorems for these quantities. Then, in Section |3.2| we show how to obtain the threshold for the 
construction of CoPE sets by an efficient multiplier bootstrap. We compare it with a method for Gaussian 
smooth noise based on Taylor and Worsley (2007). Section 3.3 combines the previous results in a concise 
algorithm for the construction of CoPE sets. 

Section [4] shows a toy example to investigate the non-asymptotic performance of CoPE sets and finally, 
in Section [5] we apply our methods for a general linear model to the climate data. All proofs are in the 
appendix. 


2 Error control for excursion sets - CoPE sets 

The domain S C on which all our functions and processes are defined, is assumed to be a compact but 
not necessarily connected subset of Euclidean space. We call the topological boundary dA c of the excursion 
set A c the contour of g, at the level c. 

Assumption 1. We assume that 

(a) the target function /i is continuous and the level set {s : /x(s) = c} is equal to dA c . 

(b) the estimator jx n { s) is continuous in s (for all n £ Nj. 

(c) there is a sequence of numbers r„ and a continuous function a : S —> R + such that 

7 -> goo, 


( 2 ) 
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weakly in C(S). Here, G is a Gaussian field on S with mean zero, unit variance and continuous sample 
paths with probability one. 


We will obtain nested estimates by thresholding the surface fi n { s) as follows: 

fin{s) - c 


A+ := ic(+a) := <j s : ° > +a 

T„cr(s) 


A c : A c ( a) 


r(s) 


> —a 


( 3 ) 


where a + and a - are appropriate non-negative constants to be determined. Note that in this notation 
A c = A c (0) and A^ are themselves excursion sets. Moreover, for any choice of a > 0 we have the inclusions 
Afi C A c C A~, and hence the estimates obtained via (p| are in fact nested. The function used to define 


the excursion sets is similar to the test statistic used in French and Sain (2013). 

The following main result shows how the constant a in ([3j) can be chosen such that At Ci c C A~ with 
a predefined probability. 


Theorem 1. If the Assumptions^ hold, then 

[At 


lim P 

n—> oo 


c A c c A: 


= P 


sup |G(s)| < a 
dA c 


A direct consequence of Theorem [l] is that At C A c C A~ with asymptotic probability at least 1 — a if 
we choose a such that P [sup sgaj 4 c |G(s)| > a] < a. The determination of a poses a computational challenge 
since the distribution of the supremum of |G(s)| and the set dA c are unknown. In Section 3.2 we propose 
an easy and fast way to approximate this distribution by a multiplier bootstrap. 

As mentioned in the Introduction, confidence sets for the excursion set A c yield confidence sets for the 
contour dA c . More precisely, we have the following 

Corollary 1. Under the assumptions of Theorem [7| we have 


lim P 

n—> oo 


dA c c cl ( A 


(A- \ i+) 


= p 


sup |G(s)| < a 
,dA a 


where cl denotes the topological closure. 




Note that, conversely, confidence sets for the contour do not automatically give confidence sets for the 
excursion set. In Figure [ 2 ] we show a simple schematic example of a pair of nested sets Afif for which 
dA c C Af \ A c holds but Af C A c C A~ does not. 

In fact, excluding these cases is the more labori¬ 
ous part of the proof of Theorem [T| The key is to di¬ 
vide the region S into a close-range zone where /i(s) 
is close to c and a long-range zone. More precisely, 
the close-range zone is given by the inflated bound¬ 
ary A^ = {s € S : c — i]cr( s) < p(s) < c + r]a(s)}. 

Then, the strategy of the proof is to let the param¬ 
eter 77 go to zero at an appropriate rate as n -) 00 
such that, eventually, the probability of a part of 
Af falsely appearing in the long-range zone S\A% 

(as shown in Figure [2| vanishes. The probability of 
making an error remains in the close-range zone, and 
is asymptotically given by P [sup aj4c |G(s)| > a]. 

We want to emphasize that Theorem 1 and its 
Corollary are valid for any estimator fi n satisfying 
Assumption 1. Thus, they hold generally whether 

the estimator is based on an increasing number of repeated observations (like in our data example) or an 
increasing number of sampling spatial points (like in the spatial statistics and nonparametric regression 
problems). However, for concreteness, we focus on the former situation, which we develop in detail in the 
following section. 


Fig. 2: A simple example of nested sets At that bound 
the contour dA c but do not satisfy the inclusion 

At C A c C At 
























3 CoPE sets for general linear models 


5 


3 CoPE sets for general linear models 

For concreteness and application to the climate data, we here present how CoPE sets are obtained, in theory 
and in practice, when the target function is a parameter function in a general linear model and n is the 
number of repeated observations. 


3.1 Asymptotic coverage probabilities 

We begin by proving an analog of Theorem[l]for the parameters of a general linear model. The most difficult 
part is to establish a Central Limit Theorem as in ([2]). This will require conditions on the error field as well 
as on the design. We consider the model 

Y (s) = Xb(s) + e(s), seSd 1 (4) 


where Y(s) is an n x 1 vector of observations, X is a known n x p design matrix, b(s) = (&i(s) 
is an unknown px 1 vector of parameters and e(s) = (ei(s),..., e n (s)) with e±,. . . 


>M S )) 

an unknown 


stochastic process. Results of the kind presented here are well-known (see e.g. Eicker (1963)). We show and 
prove versions tailored for our specific purpose for coherence and convenience. 

The least squares regression estimator for b(s) in the model Q is 


-l. 


b(s) = (X 1 X) X 1 Y(s). 


In the notation of Section ([2]), the target function p is now one of the parameter functions of the model Q, 
b i, say. Of course, the choice of b\ is arbitrary and any other coefficient of b may be considered, with the 
obvious modifications of the assumptions and theorems. Naturally, 6i(s) now plays the role of the estimator 
An(s). 

Further, define a : S — > R> 0 via <r 2 (s) = var [e(s)] and the correlation function c : S x S —> (—1,1) by 


c(si,s 2 ) 


cov[e(s 1 ),e(s 2 )] 

ct(si)ct(s 2 ) 


Recall that for 1 < p < oo, the p-norm ||A|| p of a matrix A is defined to be ||A|| p = sup|| a ,|| p=1 ||Aa;|| p . 
Hence, by definition ||Ax|| p < ||A|| p ||a;||p for all x. In the special case p = oo the matrix norm II^Hoo is the 
maximum absolute row sum of A , i.e. 


11^4- — (^q )l<i<m ,l<j<n | |oo — max \ ' dj.j | ■ 

Ki<m z ' 

- “ 3 =1 

Definition 1. (a) For vectors s,t £ R N define the block (s,t] = x ••• x (snAn] C and for 

a stochastic process e(s) with index set containing (s, t] define the increment of e(s) around (s, t] (cf. 
Bickel and Wichura] ([l97i|) as 

e (( s > t D= e( Sl + Ki(ti - Si), •• • ,s N + K N (t N - Sat)) . 

—0,1 K, iV —0,1 


(b) We denote the Lebesgue measure of a set A C S by |H|. For non-negative numbers 6, 7 , /? we say that 
the error field e(s) has the properties 

• Nl-<5, if sup sgS cr(s) - ( 2 + ' 5 ).E|e(s )| 2+,5 < 00 ; 

• N 2 -( 7 ,/ 3 ), if there exists a constant C > 0 such that £ , |(tT _ 1 e)(R )| 2+7 < C\B\ 1+l3 for all blocks 
B c S. 

In dimension N = 1 the definition of an increment yields e((si,ti] = e(ti) — e(si)) the usual increment. 
In Dimension N = 2 we get e((s, t]) = e(ti, t 2 ) — e(si, f 2 ) — e(ti, s 2 ) + e(si, s 2 ). 


Assumption 2. Assume that 
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(a) the parameter functions b(s) = (&i(s),..., b p (s)) are continuous and the level set {s : &i(s) = c} is 
equal to dA c (bi). 

(b) the noise field e(s) has continuous sample paths with probability one. Moreover, a centered unit variance 
Gaussian field with correlation function c(si,S 2 ) also has continuous sample paths with probability one. 


(c) the variance function <r(s) is continuous. 

(d) there exists a 5 > 0 such that e(s) has the property Nl-<5 and n | |X(X T X) _1 / 2 | —> 0 as n —> oo. 

(e) there exist 7 > 0 and fi > 0 such that e(s) has the property N 2 -( 7 , /3) and max n6 pj n | |X(X T X) _1//2 | |^~ 7 < 
00. 


Part (a) and (b) of Assumption [ 2 ] are tantamount to the first two conditions in Assumption [l] Parts 
(c), (d) and (e) will ensure that the parameter function bi enjoys a Central Limit Theorem. Note that the 
assumptions on the increments of the error field and on the design matrix X are coupled. The following 
Theorem [2] gives convergence results and explains how we can obtain CoPE sets for b±. 

Theorem 2. Under Assumption^ the following is true. 

(a) the weak convergence 

\/X T X (b(s) - b(s)) /ct(s) -> G 0p (s), 

holds, where G® p is an W-valued mean zero, unit variance Gaussian random field with correlation 
function 

cov [G® p (s 1 ),G 0p (s 2 )] = c(s!,s 2 )Ip. 


(b) if additionally the top-left entry ir n = 


[X T X 


-1 


(1,0,..., 0), the first standard basis vector, 


J11 


of the matrix X T X is not zero and with ej = 


M V 2 ef(X T X)-V 2 


v T e 


then we have 

IMI2 1 M v M s )~ 1 (>i( s ) - M s )) -> G ( s )> 

weakly, where G is a mean zero, unit variance Gaussian field on S with correlation function cov [G(si), G(s 2 )] 
c(si,s 2 ). 

(c) under the additional assumptions of part (b), and if we define 


then 


:= < s : 


&i(s) - c 

in V 2 / \ 

|v|| 2 7 if cr(s) 


> ±a 


lim P 

At(h)c A c {b 1 )cA-(b 1 ) 

= p 

sup G(s) < a 

n—± 00 



_< 9 A c (bi) 


( 5 ) 


3.2 Approximating the tail probabilities of G 

3.2.1 Multiplier bootstrap 

In order to obtain CoPE sets from Theorem [2] we need to know the tail distribution of the supremum of the 
limiting (non-stationary) Gaussian field G. In applications, as for example our climate data, the distribution 
of G (and hence of its supremum) is unknown, because it depends on the unknown (non-stationary) covariance 
function. In our motivating application the only information we have about G is contained in the residuals 
(i?i(s),... ,R n ( s)) = R(s) = Y(s) — Xb(s) of the linear regression. 

A way of approximating the distribution of the limiting Gaussian field G in this situation is given by 


the multiplier or wild bootstrap first introduced by Wu (1986) and later studied by Mammen (1992 19931, 
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Hardle and Mammen ([1993). It is based on the following idea. Let g\,... ,g n be i.i.d standard Gaussian 
random variables independent of the data. Consider the random field 


G(s) = n z^gjRji s). 


( 6 ) 


3 = 1 


Then, conditional on the residuals {i?. ? (s)}’? 1? the field G(s) is Gaussian and has covariance 


cov 


G(si),G(s 2 ) = - ^2 R l (si)R j {s 2 )cov[g l ,g j ] 


i,j =1 




(7) 


i=i 


the sample covariance. For large n, we expect the sample covariance to resemble the true covariance. The idea 
is to take the distribution of G as an approximation of the distribution of G. In particular, we can approximate 


P [sup sea4c |G(s)| < a], needed in Theorem 


B by 


p 


sup se0Ac |G(s)| < a I {Rj( s)}" 


latter can be efficiently computed by generating a large number M of i.i.d. copies Gi(s), 


In practice, the 
,,G m (s) of G(s), 


conditional on {i?j(s)}” =1 , via ([6]) and evaluating M 1 


su P s eS 


Gj( s) 


< a 


Theoretical considerations: The above approximation of the distribution of the supremum requires 
some justification because the sample covariance ([7]) itself is not a good estimator of the true covariance 
function in our high dimensional setting, where the number of locations is much higher than the sample size 
n (about ten thousand grid points vs. 58 field realizations in the climate data). However, the claim is about 
the distribution of the supremum of the process instead. For a discrete set of locations, the distribution of 
the supremum sup^g^ |G(s)| is similar to the distribution of the maximum of a high-dimensional Gaussian 
random vector, recently considered by Chernozhukov et al. (20131. Substantially extending the results of 


Mammen (19931 in the high dimensional setting, Chernozhukov et al. (20131 show that the distribution of 


the maximum can be well approximated by the Gaussian multiplier bootstrap using realizations of a not 
necessarily Gaussian random vector with the same covariance matrix. In this sense, the multiplier bootstrap 
is valid in our setting. This is confirmed by simulations in Section 4 below. 

Computational considerations: Besides these theoretical considerations, the multiplier bootstrap is 
also computationally attractive. In comparison with the direct simulation of the limiting field. While it is 
theoretically possible to simulate a number of realizations of a non-stationary Gaussian field with a given 
covariance as in Q and to obtain tail probabilities from these, in practice this is computationally infeasible. 
Assuming that all fields are observed at L locations in S the direct method first requires computing the 
L x L covariance matrix which is of complexity 0(nL 2 ). Then, a Cholesky decomposition of the covariance 
matrix must be computed in time 0(L 3 ). Finally, for each realization a matrix-vector product with the 
triangular matrix from the Cholesky decomposition is required and hence N realizations can be obtained in 
0(L 2 N) time. This yields a total complexity of 0(L 3 + L 2 ( n + N)) which is prohibitively large (in our data 


L = 9051). This problem has also been encountered by Adler et al. (2012) in a setting where information 
about the field is available through the true covariance function instead of realizations of it. 

In contrast, creating one multiplier bootstrap realizations of the field G requires computing a linear 
combination of n vectors of dimension L which is of complexity 0(nL). Hence, N multiplier bootstrap 
realizations can be generated in O(nNL) time. Further, note that the simulation of N bootstrap realizations 
can be written as a matrix multiplication. Let E be a L x n matrix such that each column is one residual 
Rj and let V be a n x N matrix with i.i.d. standard Gaussian entries. Then, the columns of EV correspond 
to realizations of G. This makes the multiplier bootstrap very efficient because very fast implementations of 
matrix multiplication are available and it is an operation that can easily be parallelized 

A comparison of the computational burden to existing methods is not easy, since all work that is known to 
the authors considers different settings. However, to put the above considerations in perspective, we briefly 


discuss computational requirements of the method proposed by French (2014). This method requires (French 
2014 Sec. 3.7) the computation of a kriging estimate in time 0{n 3 ), where ni is the number of observed 


locations, and a Cholesky decomposition which is of complexity 0((ni + n g ) 3 ), where n g is the number of 
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grid locations (corresponding to L in our notation). This increased complexity compared to the multiplier 
bootstrap is reflected in actual empirical computation times. According to his own experiments, the method 


of French (2014) applied to n ~ 100 observed locations and a grid size of 100 x 100, can be computed, on 


average, in just under five minutes. In contrast, the entire data analysis for the mean summer temperature, 
including linear regression at each grid point and computation of CoPE sets with the multiplier bootstrap, 
with a grid size of L ss 10, 000 and n = 58 is performed in under five seconds on a machine comparable to 


the one used by French (20141. 


3.2.2 An alternative method for smooth fields 

If one is willing to assume that the limiting field is twice differentiable and that the error field e(s) in p} is 
Gaussian itself then the Gaussian Kinematic Formula (GKF) (e.g. Taylor (20061, Taylor et al. (20051, Taylor 


and Worsley (2007)) offers another way to approximate tail probabilities. Our motivation for presenting it 
here is twofold: First, it offers, under the additional assumptions made above, an elegant and accurate 
way of computing tail probabilities of Gaussian fields that does not require simulations and that has been 
successfully applied (Taylor and Worsley 2007). Second, it shows that, at least for smooth fields, the tail 
probability of the supremum is intrinsically low-dimensional. More precisely, for a field on it is given 
(up to an exponential error term) by N +1 numbers. This gives an additional justification for the ability of 
the multiplier bootstrap method to estimate the tail probability despite the high dimensionality of the field. 

The GKF is based on two properties of smooth Gaussian fields. The first is that for such fields the 
exceedance probability for high thresholds can be approximated by the expected Euler characteristic. More 
precisely, for any set B C S with smooth boundary, 


P 


sup G(s) > a 

s gb 


= E \x(A a (G) n B)\ + G(exp(-a 2 /2)), 


( 8 ) 


where \ is the Euler characteristic (see e.g. Adler and Taylor (2007)). The second property is that the 
expected Euler characteristic in ([8]) has a closed formula. Indeed, with A(s) = var G(s) , where G(s) is the 
vector of partial derivatives of G, we can write (cf. e.g. Taylor (2006)) 

N 

E [ X (A a (G) nB)}=J2 £d{B, A)pd(a), (9) 

d —0 


with Cd the d-th order Lipschitz-Killing curvature (LKC) (see e.g. Taylor (2006) for details on these quanti¬ 
ties) and known functions pd{a). We can use this to obtain CoPE sets: If the Assumptions ([!]) are satisfied 
then Theorem ([!]) implies in conjunction with ([8]) and|9]) that 


lim P 

n—too 


At c A c c A„ 


= 1 -p 

N-\ 


sup |G(s) | > a 

s £dA c 


= 1 - Cd(dA c , A )p d {a) + £>(exp(-a 2 /2)). 


( 10 ) 


d=0 


To use (101, the problem amounts to estimating the LKCs of the boundary dA c . Taylor and Worsley (2007) 


propose a method to estimate the LKCs based on a finite number of realizations of G. Applying this method 
requires a triangulation of the plug-in estimate dA c of the boundary dA c . 


While the triangulation is challenging, yet feasible, Taylor and Worsley (2007) prove the validity of their 


method only when applied to realizations of the Gaussian field G. In our application, however, we only 
have realizations of a generally non-Gaussian field (the residuals in of the linear model, cf. Section [5]) with 
asymptotically the same covariance as G. For completeness, we compare this method to the multiplier 
bootstrap method in the simulations section below. 


3.3 Algorithm 

Combining the results of the previous sections we can give the exact procedure for obtaining CoPE sets for 
the parameters of the linear model Y(s) = Xb(s) + e(s). 




















































4 Simulations 


Algorithm 1. Given a design matrix X and observations Y(s) following the linear model 0- If Assump- 
tions^hold, the following yields CoPE sets for b\(s). 


(a) 


Compute the estimate b(s) = (X T X) 
Xb(s). With the empirical variance cr 2 (s) 
<x(s) _1 R(s). 


X t Y(s) and the corresponding residuals R(s) = Y(s) — 
= n~ l -Sy ( s ) compute the normalized residuals R(s) = 


(b) Determine a such that approximately P |sup„ p <; |G(s)| > a] < a. For example, use the multiplier boot¬ 


strap procedure presented in Section 


3.2.1 


with the residuals 


mi 


to generate i.i.d. copies of a 


Gaussian G field with covariance structure given by the sample covariance. With these, determine a 


such that P 


SU PsGS 


G(s) 


> a 


< a. 


(c) Obtain the nested CoPE sets defined in equation © 


4 Simulations 

This section includes some artificial simulations to show that the proposed methods provide approximately 
the right coverage in practical non-asymptotic situations with non-smooth, non-stationary and non-Gaussian 
noise. We will describe ways to obtain error fields with these properties. Our objective in the design of the 
error fields described below is not to imitate the data but to introduce non-stationarity and non-Gaussianity 
in a transparent and reproducible way, showing the full potential of the method. In fact, the error field 
that we encounter in the data (cf. Section ©) is better behaved as far as smoothness and stationarity are 
concerned than the artificial fields we investigate here. 


4.1 Setup 

As a simple instance of the general linear model Q we consider the signal plus noise model 

%(s) =/x(s) + ej(s), j = 1,... ,n 

with non-stationary and non-Gaussian noise e(s) over a square region S of size 10 x 10 consisting of 64 x 64 
square pixels. We will consider three different noise fields and we describe in the following how to obtain a 
realization of each. 


Noise 1 In the upper half of S, each pixel is assigned the value of a standard normal random variable, all of 
which are independent. In the lower half, the pixels are grouped together in blocks of 4 by 4 pixels and 


each block is assigned the value of a standard normal, again all independent (cf. Figure 3b). Finally, 
the entire picture is convolved with a Gaussian kernel with bandwidth one and all values are multiplied 
by a scaling factor of 50. 

Noise 2 Identical to Noise 1 except the image is smoothed by a Laplace kernel with bandwidth one instead 
of a Gaussian and the scaling factor is 100. 


Noise 3 Each pixel in the upper half is assigned the value of a Laplace distributed random variable with 
mean zero and variance two. In the lower half, pixels are assigned the values of independent Student 
t-distributed random variables with 10 degrees of freedom. The entire picture is convolved with a 
Gaussian kernel of bandwidth one and multiplied by a scaling factor of 25. 


The noise fields Noise 1-3 are intentionally designed to have non-homogeneous variance and scaling factors 
are chosen ad-hoc such that all three fields can be conveniently displayed on a common scale. 

The signal p is a linear combination of three Gaussians. Figure[3a]shows the signal p. Each one realization 
of the three noise fields is shown in Figure [3] 

We controlled the probability of coverage at the level 1 — a = 0.9 using Theorem |T] The estimator for p 
here is the mean p n { s) = n _1 ^' l =1 Vj( s ) an d the thresholds for the CoPE sets are obtained using Algorithm 

[llwith 11 v 1 1 1 7Tn ^ = \Jn. The threshold a was computed using the multiplier bootstrap procedure proposed 
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(a) The signal fi(s) with 
the true contour 
dA c (fi) for c = 4/3 in 
purple. 



(b) A realization of Noise 
1 before smoothing 
with a Gaussian and 
scaling. 


I 



0 2 4 6 8 10 


(c) Noise 1. Nonstation¬ 
ary Gaussian noise 
smoothed with a 
Gaussian kernel. 



(d) Noise 2. Nonstation¬ 
ary Gaussian noise 
smoothed with a 
Laplace kernel. 



0 2 4 6 8 10 


(e) Noise 3. Laplace and 
Student-t distributed 
noise smoothed with a 
Gaussian kernel. 


Fig. 3: The signal and noise in the toy example. 


using either the true boundary dA c or the plug-in estimate dA c . The results of our method 


in Section [3.2.1 

using dA c for each one run with the three noise fields and sample sizes n = 60, n = 120 and n = 240 are 
shown in Figure [4j 


4.2 Performance of CoPE sets 

We analyzed the performance of our method on each 5000 runs of the toy examples shown in Section |4.1| 
with sample sizes n = 60, n = 120 and n = 240. Table [T] shows the percentage of trials in which coverage 
A+ C A c C A~ was achieved, if either the true boundary dA c or the plug-in estimator dA c was used to 
determine the threshold. 

We see that the empirical coverage is smaller than the nominal level in all experiments but approaches the 
nominal level reasonably fast as the sample size increases. In fact, when n = 240, the simulation confidence 
interval cover the nominal level of 90%, suggesting asymptotic unbiasedness. Comparing the two columns, 
we see that the non-asymptotic bias is not caused by the lack of knowledge of the true boundary. It may be 
a consequence of the bootstrap procedure instead. 


Computational performance As already noted in Section 3.2.1[ the multiplier bootstrap allows for a very 
fast computation of CoPE sets. In the simulations, the CoPE sets for a sample of size n = 240, each on a 
grid of 64 x 64 = 4096 locations could be computed in less than two seconds on a standard laptop. 


4.3 Comparison with Taylor’s Method 

In this Section we compare the multiplier bootstrap with the method proposed by Taylor and Worsley (2007), 
as described in Section 3.2.2 We use both methods to approximate the distribution of sup sgaj 4 c |e(s)/cr(s)|, 
where e(s) is distributed according to Noise 1 (see Section 4.1 above), cr 2 (s) = var[e(s)] and dA c is the 


contour A c (fi',) of the function /x shown in Figure |3a| at level c = 4/3. The true cumulative density function 
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n = 60 


n = 120 


n = 240 



Fig. 4: The output of our method for the three noise fields described above (corresponding to rows) and 
for sample sizes n = 60,120, 240 (corresponding to columns) with the target function /i(s) shown in 
Figure 3a In all pictures we show a heat map of the estimator /t n (s), the boundary of A c (/x) in 
purple as well as the boundaries of %+ and A~ in red and green, respectively. The threshold a was 
obtained according to Theorem[l]to guarantee inclusion A+ C A c C A~ with confidence 1 — a = 0.9. 



dA c 

dA c 

Noise field 1 

n = 60 

86.20% ± 0.49% 

86.16% ±0.49% 

120 

88.62% ± 0.45% 

88.74% ± 0.45% 

240 

88.94% ± 0.44% 

88.9% ± 0.44% 

Noise field 2 

n = 60 

87.22% ± 0.47% 

88.74% ± 0.45% 

120 

89.22% ± 0.44% 

89.26% ± 0.44% 

240 

89.76% ± 0.43% 

89.70% ± 0.43% 

Noise field 3 

n = 60 

86.44% ± 0.48% 

86.62% ± 0.48% 

120 

88.60% ± 0.44% 

88.78% ± 0.45% 

240 

89.76% ± 0.43% 

89.94% ± 0.43% 


Tab. 1: Percentage of trials in which A+ C A c C A c . The nominal coverage probability is 90%. 
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- True distribution 
-Multiplier Bootstrap 
Taylor's Method 


Fig. 5: The probability P [supggg^ |e(s)/cr(s)| < a] (the horizontal axis shows a) and approximations of it 
via the multiplier bootstrap or Taylor’s method based on a sample of size n. Here, the error held e(s) 


has the distribution described in Noise 1 and A c = A c (n) for the function /r shown in Figure 3a 


for supgg^^ |e(s)/cr(s)| and its empirical approximations based on the multiplier bootstrap and Taylor’s 
method are shown in Figure [5] The empirical cdfs are each based on a single i.i.d. sample ei(s),... ,e n (s) 
for n = 10, 30 and 60. For the multiplier bootstrap we generated 5, 000 bootstrap realizations. The true cdf 
was calculated empirically using 10,000 i.i.d. samples of e(s). 

Both methods give a remarkably good approximation of the true distribution of the supremum, particu¬ 
larly for sample sizes of n = 30 and higher. However, while Taylor’s method only gives a valid approximation 
in the tail of the distribution, the multiplier bootstrap approximates all parts of the cdf. 

5 Application to the climate data 
5.1 Data setup 

In our application we have a total of n = n O) + n^ observations, the first nO) observations are the ’past’, 
the last n < - h> are the ’future’. Within each period we model the change in mean temperature linearly in time. 
More precisely we have 


Y,-(s) = T (q) (s) + m (o) ( s)^ o) + €j (s), j = • ,n (a) 

Yj(s) = T < - b \ s) + m^(s)tf' > + ej( s), j = + 1,... ,n^ + \ 




( 11 ) 


Without loss of generality, we may assume that J^" =1 tj = 0 an d X)J= n (“” + i tj = 0. We will denote 
the covariance of the error field e(s) by c(si,S 2 ) = cov [e(si), e(s 2 )]. Our goal is to give CoPE sets for the 
excursion sets of the difference T^ b \ s) — T^ a \ s). Therefore, we define the parameter vector and the design 
matrix 


/ \ 

( 

b 2 (s) 


h(s) 


\ h(s) ^ 

V 


T(°)( s) 
m.(°) (s) 

(s) 


in 


( b) 


X = 


( o 

0 

1 


1 t 


(a) 


\ 


C nC“> 

0 


1 1 


0 


Ab) 

(, ra (“)+1 


f (b) 


to be able to rewrite (111 as a general linear model Y(s) = Xb(s) + e(s). Our objective is now to formulate 
Assumptions on the design and the noise under which we can apply Algorithm [T] to the data. This is done 
in the following. 
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Assumption 3. Assume that 


(a) the parameter functions b are continuous and and the level set {s : &i(s) = c} is equal to dA c (b\). 

(b) the noise field e(s) has continuous sample paths with probability one and a centered unit variance 
Gaussian field with correlation function c(si,S2) also has continuous sample paths with probability one. 

(c) the variance function er(s) is continuous. 

(d) there exist numbers 5, (3 > 0 and 7 > 0 such that the error field e(s) has the properties N1 and 

N2-( 7 , f3). 

(e) n = n^ = n /2 and that both sets of design points t^ and t^ are equally spaced (possibly with 
different spacing for the periods (a) and (b)). 


The next final and central statement now asserts that these Assumptions are indeed sufficient for Algo¬ 
rithm [T] to be valid. Its proof is a direct application of Theorem [2j 

Proposition 1. Under Model (111 Assumptions [5] imply that Assumptions^ hold for the target function 
b\ (s) and the estimator bi (s) with r n = 2n~ 1 / 2 . In particular, 


2^s) (b^) - bi(s)) ->G(s), 

weakly, where G is a mean zero unit variance Gaussian field with correlation function cov [G(si), G(s 2 )] = 
c(si, S2). Consequently, Algorithm^ can be used to obtain CoPE sets for the excursion sets ofbi(s). 


5.2 Data analysis 

The results for the climate data described in the Introduction, shown in Figure [I] correspond to CoPE sets 
for 61 (s) = T^(s) — T( a )( s) obtained via Algorithm [l] with 11v11 27 r^/ 2 cr(s ) —1 = 2 /y/n and n = 29 + 29 = 58. 
The target level is c = 2 °C and nominal coverage probability is fixed at 1 — a = 0.9. 

For the mean summer temperature, it may be stated with 90% confidence that the Rocky Mountains 
and the Sierra Madre Occidental mountains of Mexico are at risk of exhibiting a warming of 2 °C or more in 
the given time period, while the Florida Peninsula, parts of the Mexican Gulf, large parts of the Canadian 
Northwest and the northern part of the Labrador Peninsula are not at risk. 

For the mean winter temperature, some regions around the Hudson Bay and in the Canadian Shield are 
identified to be at a high risk while a comparatively small region north of the Mexican Gulf is considered 
not at risk for extreme warming. 

For the computation time we remark that the entire analysis of one season, including the pointwise linear 
regression and the multiplier bootstrap to obtain the CoPE sets was performed in under five seconds on a 
regular laptop. 
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A Proofs 

Proof of Theorem [7] We start by showing that 


lim inf P 

n—t oo 


At c A c c A: 


> P 


sup |G(s)| < a 
dA c 


( 12 ) 


For r] > 0 define the inflated boundary A^ = {s £ S : c— rja(s) < fj( s) < c + 7y(r(s)}. The idea of the proof is 
that, loosely speaking, points outside of become irrelevant in the limit n — > oo since their values are far 
from c and, if we let g go to zero at an appropriate rate, we finally end up with the boundary dA c . More 
precisely, we note that 


Ms)-p(s) 


> —a for all s £ A c fl A 1 } and 


£n(s) - M s ) 


T„cr(s) C T„cr(s) 

implies that /t n (s) > c — r„cr(s)a for all s £ A c and hence A c C A~. Similarly, 


> — TfT n 1 — a for all s £ A c \ A]!, 


An(s) ^ MS) 
T„a(s) 


< a for all s £ (S \ A c ) n A v c and 


//„(s) -p(s) 
r n a(s) 


< rjT n 1 + a for all s G (5 \ A c ) \ A v c 


implies Af C A c . Combining these observations, we see that Af C A c C A c holds, provided that 
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positive numbers such that rj n — > 0 and rj n T n 1 —> oo. We can then write 


p 

A+ c A c c A c 

> P 

sup 




s eA? n 


£n(s) - Ms) 


T„cr(s) 


< a and sup 

S eS\>2" 


An(s) “Ms) 


r„cr(s) 


< a + 


> P 


sup 

sGA ? 71 


An(s) -/U(S) 


r(s) 


< a 

+ p 

sup 



ses\A? n 


Ms)-n(s) 


r(s) 


< a + r? n r n 


- 1 . 


(13) 


(P 


UP 


We first show that the term (II) goes to zero. To this end let <5 > 0 arbitrary. Let b £ 1 such that 
P [sup s6S |G(s)| < b] > 1 — 6 and no £ N such that a + > b for all n > no- Also, let n± large enough 

such that 


p 

sup 


,sG5 


M s) ~ p(s) 
t„ct(s) 


< b 


- P 


sup |G(s)| < b 
_s es 


<S, 


for all n>n\. In consequence, for all n > max{?ro,ni} 


sup 

M s ) - M) 

< a + rj n T n 1 

> P 

sup 

M s ) - M) 

< b 

t„<t(s) 

r„cr(s) 

ses\A 2" 



.seS 



>1-2 5 . 


Since 6 > 0 was arbitrary it follows that (II) converges to zero as n —> oo. 

To prove convergence of (I) we need the following 

Lemma 1. Under Assumptions^ part (a) if rj n —> 0 then the Hausdorff distance 5 n := dn (Aj?", dA c ) —> 0. 

Proof. Let us define the set ( dA c ) e := {s £ S : d(s,dA c ) < e}. We prove the assertion by showing that for 
any e > 0 there exists an ?y > 0 such that A% C ( dA c ) e . To this end, assume the contrary. Then, there 
exists e > 0 such that for any > 0 we find s v £ A^ with d(s r ,,dA c ) > e. The sequence (s^p) is contained 
in the compact set S and hence has a convergent subsequence with limit s*, say. By construction, we have 
s* £ n r; >o A]! = dA c . On the other hand, 0 = d(s*,dA c ) = lim^o d( s, ; , dA c ) > e, a contradiction. □ 

Recall that for a function / : S —> R. and some number <5 > 0 the modulus of continuity is defined as 
w{f,S) = sup| Sl _ S2 |< 5 |/(si) - f (s 2 )|. Since {^^(s)- 1 (£„(s) - M))} neN is weakly convergent, we have 


lim lim sup P 

ft —*^0 n—^oo 


W 


An(s) ~ /i(s) 
T„a(s) 


,6 



= 0 


(14) 


for all positive £ (Khoshnevisan 2002 Prop. 2.4.1 and Exc. 3.3.1). Together with Lemma[l]this implies 
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in probability, 
yields 


Since sup sgaAc |(/i„(s) — /i(s))/(r ra cr(s))| converges in distribution to sup se9j4r |G(s)| this 

M s ) - /x(s) 


sup 

seA2 n 


r(s) 


sup |G(s)| 
s EdA c 


in distribution. In view of (131 this completes the proof of (12). 
It remains to prove the opposite inequality, i.e. 


lim sup P 

n—> oo 


A+ c A c c A c 


< P 


sup |G(s)| < a 
_dA c 


(15) 


If for some arbitrary S > 0 we have r n 1 cr(s) 1 (//„(s) — c) > a + S for some s £ dA c then by continuity there 
is a s £ S \ A c for which t“ 1 <t(s ) _1 (/z n (s) — c) > a and hence the inclusion Af C A c does not hold. Since 
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an analogous argument works for the inclusion A c C A c , we have 


At c A c c A„ 
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Since S > 0 was arbitrary and supgg^ |G(s)| has a continuous distribution the bound (151 follows. □ 


Proof of Corollary [7] For any pair of nested sets Af we have that Af C A c C A~ implies dA c C cl(A“ \ Af). 
On the other hand, the latter will certainly fail to hold if sup sgayic 
observations yields 


Ar»(s)~C 
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At C A c c A r 


< P 


dA c cd(A-\At) 


< P 


> a. Combining these two 
An(s) ~ C <a 

□ 


sup 
sEdA c 


Taking the limit n —> oo of this inequality and using Theorem [T] gives the assertion. 

Proof of Theorem [1[ We begin by proving part (a). Let us define A = (X T X) 1 X T , giving b(s) — b(s) = 
Ae(s). In order to prove weak convergence of the process we first show convergence of the finite dimensional 
distributions and then tightness of the sequence (Khoshnevisan 2002) Prop. 3.3.1). 

For the former, let si,..., sk £ S be arbitrary. We need to show that with 

U = (TFX^Si)- 1 (b( Sl ) - b( Sl )) ,..., v / X^X<t(sk) _1 (b{s K ) - b(s*))) 


= I 


Vx T x 


(cr(si) 1 Ae(si),..., cr(sic) 1 Ae(s^)), 

(here, '®' denotes the Kronecker product of two matrices) we have convergence U —> A/”(0, [c(sj, Sj)]fj = 1 < 8 >I p ) 
in distribution. We readily see that E[U] = 0 and for the covariance we compute 

cov [U] = I K ® VX T XA E | (cr(si) _ 1 e(si),..., a(s K )~ 1 e{s K )) (er(si) _ 1 e(si),..., er(sK) _ 1 e(s;f)) T } 


x I 


a t Vxjx 

= {lie ® [Vx^xa] } |[c( Si , s,)]5 =1 <g> I„} [l K <g> [aVx^x] } 
= [c(8i, s,)]5 =1 ® Vt^XAA T Vt^X = [c(si, Si )]5 =1 ® I p . 


v L v L1 3 = 

We employ the Cramer-Wold device to show convergence of U. Indeed, let and (aq,..., ape) £ be 

some fixed arbitrary vector and compute 


(U,a) = ((a{ si) 1 e(si),...,o-(s if ) ^(s^)) ,1k ® [a 7 VX t X a) 

K K n 

= ^a(s i )- 1 (e(s i ), AVx^Xa;) = a(s l )- 1 e,(s l ) ^A T VX T Xa, 

i= 1 z=l ji=l 

By interchanging the sums and defining Wj = ttiLi a ( s i)~ le ji s i) ( 'A T VX T Xai ^ , we have managed to write 
(U, a) as a sum of independent random variables (U, a) = JT Wj. The goal is now to use the CLT in the form 
of Lyapunov for the random variables Wj. To this end compute var t2’j=i Wj = var [(U, a)] = a T cov [£/] a 
and note that since we have already showed U to have the right covariance the claimed convergence will 
follow once we establish the Lyapunov condition. For this purpose let S be as in Assumption [2] to give 

K 2+5 


E W +5 < Si ) 

3= 1 3 =1 *=1 

n K 


j=i»=i 


Vx T x« t ) 


| 2+<5 


2+<5 


A T VX T Xai 


< CK 2 \\a\\ no n 


A T Vt^X 


2+5 
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as n —> oo. This concludes the proof of convergence for the finite dimensional distributions. 

It remains to show tightness of the sequence | \/X T XAcr(s) _ 1 e(s)|. For any block B C S we have with 
7 and /3 as in Assumption [2] 


E 


v / X T XA( f 7 - 1 e)(B) 


2+7 


< 


2+7 

VX T XA “ 5Zs|(<r _ 1 e)(B)p 


< Cn 


j=i 

2+7 


A T VX T X ' X(B) 1+ P < C'X(B) 1+0 , 


which implies tightness (cf. Bickel and Wichura (1971 Thm. 3)). 

The statement (b) follows from the fact that 

7 r n V2f7 ( s ) _1 (M s ) ^ M s )) = Kn 1/2e i {X T X)- 1/2 a{s)~ 1 VX T X (b(s) - b(s)) 

—f v t G® p (s) = ||u|| 2 G(s), 

where the weak convergence from the previous part is used. 

The last part (c) of the Theorem is a direct application of Theorem [lj where parts (a) and (b) guarantee 
that the assumptions are satisfied. 

Proof of Proposition [7} In order to be able to apply Theorem [2] we compute 


□ 


X T X = 


/ n / 2 n / 2 0 0 \ 

n /i n 0 0 

0 0 w (a) 0 

\ 0 0 0 oj W 




,(“) — 


e (<r). 


n‘“>+n< b > 


3 =1 


E (*«)■ 


,(b) 


j=n<“) +1 


It follows that 


(X T X)~ = 


( 4/n - 2 A 


0 \ 
0 


0 

— 2 /n 2/n 0 

0 0 i/“ (a) 0 

V o o o V“ (a) 7 


, ( xTx )- 1/2 = Vb 


( 6/n —2/ n 

—2/n ^/n 

0 0 


V 


0 0 


10 

ntj( a ) 

0 


o \ 

0 

0 

nuj ( ~ b '> / 


With this we obtain 


l|x ( x T x)—< A + 

v« vuA“) 

Now we note that since the design points t[p and t^ are equally spaced by Assumption ji] we have 
maxi<j< n \t^ | = 0(n) and u/ a ) = 0(n 3 ), and the same is true for the (^-counterparts. This shows 
that ||X(X 2 X)" 1 / 2 !!^ = 0(n -1 / 2 ) and therefore Assumptions [ 3 ] imply Assumptions [ 2 J Now, in the notation 
of Theorem [ 2 ] we have 7 r„ = V 7/2 and 

yfn pn 


1 /*epX T X)- 1 / 2 = ^lJ-( 3 (6 — 2 ) =: v T , 

1 \ > 9 A/ IQ n o,/Tn ' > ’ 


2\/l0 


so that 11v|| 2 = 1. This finally gives ^cr(s) 1 ( 6 1 (s) — 61 (s)) -+ G( s). 


□ 






































